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ABSTRACT 

We provide a new framework for the joint analysis of cluster observations (JACO) using simultaneous 
fits to X-ray, Sunyaev-Zel'dovich (SZ), and weak lensing data. Our method fits the mass models 
simultaneously to all data, provides explicit separation of the gaseous, dark, and stellar components, 
and — for the first time — allows joint constraints on all measurable physical parameters. JACO includes 
additional improvements to previous X-ray techniques, such as the treatment of the cluster termination 
shock and explicit inclusion of the BCG's stellar mass profile. An application of JACO to the rich 
galaxy cluster Abell 478 shows excellent agreement among the X-ray, lensing, and SZ data. We find 
that Abell 478 is consistent with a cuspy dark matter profile with inner slope n = 1. Accounting for 
the stellar mass profile of the BCG allows us to rule out inner dark matter slopes n > 1.1 at the 99% 
confidence level. At large radii, an r~ 3 asymptotic slope is preferred over an r~ 4 behavior. All single 
power law dark matter models are ruled out at greater than the 99% confidence level. JACO shows 
that self-consistent modeling of multiwavelength data can provide powerful constraints on the shape 
of the dark profile. 

Subject headings: Dark matter — X-rays: galaxies: clusters — gravitational lensing — galaxies: clusters: 
individual (Abell 478) 



1. INTRODUCTION 

N-body simulations of our Cold Dark Matter (CDM) 
dominated universe suggest that evolved galaxies and 
clusters of galaxies are contained in dark halos with 
self-similar density profiles. In the collisionless limit, 
pdm(f / fd) / S = f(x). The characteristic density S and 
radius r d depend on the halo mass and the cosmological 
parameters, but the normalized profile f(x) itself does 
not (Navarro, Frenk, & White 1997; Moore et al. 1998; 
Ghigna et al. 2000; Subramanian, Cen, & Ostriker 2000; 
Fukushige & Makino 2001; Williams, Babul, & Dalcan- 
ton 2004; Merritt et al. 2005). The robustness with 
which collisionless N-body simulations predict a univer- 
sal profile makes the observational measurement of f(x) 
particularly important. If we are to be satisfied that our 
understanding of hierarchical structure formation is cor- 
rect, we need to use the best available data to test the 
theoretical dark matter profiles in detail. 

Observations of f(x) also have implications for self- 
interacting dark matter (SIDM) as an alternative to the 
collisionless models. If the weak interaction cross sec- 
tion of the particles in a halo is large enough, substan- 
tial effects on the shape of the density profile should 
appear within a Hubble time. Initially cuspy profiles 
should evolve constant-density cores with a size ~ 3% 
of the virial radius (Spergel & Stcinhardt 2000; Burkert 
2000). It is an open question whether the SIDM cores 
are stable and long-lived enough to be observable (Bal- 



berg, Shapiro, & Inagaki 2002; Ahn & Shapiro 2005) 
or whether the cores quickly collapse, forming an r 
isothermal cusp (Kochanek & White 2000). If the cores 
are stable, measurements of f{x) can provide constraints 
on the interaction cross section and hence on the nature 
of dark matter itself. 

Rich, > 10 14 M Q clusters of galaxies are excellent lab- 
oratories for the measurement of f(x) because of the 
large number of independent techniques available for 
mass measurement. In this series of papers, we exam- 
ine the complex task of constraining f(x) simultaneously 
using the X-ray emitting intracluster medium (ICM), the 
Sunyaev-Zel'dovich decrement, and weak gravitational 
lensing maps. For the joint analysis of cluster observa- 
tions (JACO), we always begin with separate luminous 
and dark mass profiles, and calculate the observables — 
the X-ray spectrum, the SZ decrement, or the tangen- 
tial lensing shear — directly from the underlying physical 
model. In other words, our fits always take place in the 
data space, rather than in mass, density or temperature 
space. This approach allows for (1) deformation of the 
theoretical profiles to account for the appropriate instru- 
mental degradation at each wavelength and (2) the abil- 
ity to combine lensing, SZ, and X-ray data in a single 
grand total fit. 

Other works have examined the multiwavelength ap- 
proach. Miralda-Escude & Babul (1995) carried out 
the first comparison of physical X-ray and strong lens- 
ing models for a sample of three rich clusters; Squires 
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et al. (1996) extended the comparison to weak lensing. 
Zaroubi et al. (1998) and Zaroubi et al. (2001), using 
idealized X-ray and SZ observations of simulated clus- 
ters, showed that the unique axisymmetric deprojection 
of the gas density profile is possible by combining the 
two types of data. Dore et al. (2001), also using sim- 
ulated clusters, argued that the simultaneous deprojec- 
tion of SZ and weak lensing data is sensitive to the noise 
properties of the data. De Filippis et al. (2005) and 
Sereno et al. (2006) use combined X-ray and SZ data 
to place constraints on the triaxiality and inclination of 
the intracluster medium assuming that they are single- 
component, isothermal systems; they do not attempt to 
constrain the dark profile. Assuming a Navarro et al. 
(1997) mass profile, LaRoque et al. (2006) use joint X- 
ray and SZ data from 38 clusters to constrain the gas 
mass fraction. 

JACO is the first method to use X-ray, SZ, and weak 
lensing data simultaneously to constrain the shape of the 
dark matter profile in clusters. This paper describes a 
test of this method on the cluster Abell 478. Chief among 
the assumptions allowing the extraction of the dark mat- 
ter profile is that the gas in hydrostatic equilibrium with 
the overall gravitational potential. If hydrostatic equi- 
librium is strongly violated, then determination of the 
structure of the dark halo from the X-ray data alone be- 
comes difficult. Because ongoing or recent cluster merg- 
ers are common, the identification of clusters that are as 
close to prototypical (or "boring" ) as possible is valuable 
for detailed mass measurement. The difficulty of find- 
ing systems suitable for hydrostatic analysis is evident in 
the latest Chandra and XMM-Newton studies of relaxed 
clusters. 

Clusters that appeared smooth and relaxed before 
Chandra and XMM-Newton observations often contain 
features that make them unfit for equilibrium analy- 
sis. Examples include cold fronts (Vikhlinin, Markevitch, 
& Murray 2001; Bialek, Evrard, & Mohr 2002; Dupke 
& White 2003; Hallman & Markevitch 2004) or shocks 
(Markevitch & Vikhlinin 2001; Markevitch et al. 2002; 
Fabian et al. 2003; Fujita et al. 2004) in recent merg- 
ers. In clusters with X-ray substructure, the use of hy- 
drostatic equilibrium may yield incorrect results (Poole 
et al. 2006). In many otherwise seemingly relaxed clus- 
ters, heating by a central AGN is required to reconcile 
the short cooling times with the deficit of kT < 1 keV 
gas in the quantities predicted by quasihydrostatic cool- 
ing flow models (Tamura et al. 2001; Fabian et al. 2001; 
Matsushita et al. 2002; Peterson et al. 2003; Kaastra 
et al. 2004; Voit & Donahue 2005). Gas in the region 
directly affected by the AGN is unlikely to be in hydro- 
static equilibrium, and the equilibrium equations may 
yield incorrect results. 

The most recent work shows, however, that many clus- 
ters can be successfully modeled as equilibrium systems 
outside the region of influence of the central AGN. In 
these regions, the total density distribution frequently 
resembles a pure n = 1 Navarro et al. (1997) profile, 
and either a single or double /3-modcl provides a success- 
ful description of the X-ray surface brightness (Schmidt, 
Allen, & Fabian 2001; Arabadjis, Bautz, & Garmire 
2002; Allen, Schmidt, & Fabian 2002; Lewis, Buote, & 
Stocke 2003; Buote & Lewis 2004; Pointecouteau et al. 
2004; Buote, Humphrey, & Stocke 2005; Gavazzi 2005; 



Pointecouteau, Arnaud, & Pratt 2005; Vikhlinin et al. 
2006). However, in a few cases, evidence for a signif- 
icantly shallower n < 0.5 (Ettori et al. 2002; Sander- 
son, Finoguenov, & Mohr 2005; Voigt & Fabian 2006) 
exists, echoing similar results in several strong gravita- 
tional lensing studies (Tyson, Kochanski, & dell'Antonio 
1998; Sand, Treu, & Ellis 2002; Sand et al. 2004). 

The substantial progress allowed by Chandra and 
XMM-Newton data leads us to consider whether the con- 
straints on the dark matter structure could benefit from 
a refinement of the mass fitting technique and the use of 
optical and radio data in a simultaneous fit. In this pa- 
per we show that through the use of a multiwavclcngth 
modeling technique specifically aimed at constraining the 
the dark matter (as opposed to the total) mass profile, 
enhanced constraints on the slope of the dark mass dis- 
tribution in clusters are possible. In §2 we motivate the 
technique and describe it in detail. We apply the tech- 
nique to the rich cluster Abell 478 in §3, and compare 
the results with previous work and N-body models in §4. 
We conclude in §5. 

2. METHOD 
2.1. General Principles 

The purpose of JACO is to constrain the shape of 
the dark matter profile in clusters of galaxies via a sin- 
gle multi-parameter fit to X-ray, lensing, and Sunyaev- 
Zel'dovich data. The broad features of our technique are: 

Separation of the mass model into the gaseous, stellar, 
and dark components. Rather than fitting for the total 
gravitating mass, JACO splits the potential into three 
separately modeled components, thus guaranteeing their 
positivity. As a result JACO is incapable of generating 
unphysical mass profiles as is sometimes the case with pa- 
rameterized X-ray temperature profiles (Pizzolato et al. 
2003; Buote & Lewis 2004). The stellar mass contributes 
a substantial part of the gravitating mass within 3% of 
the cluster virial radius (Sand et al. 2004), an impor- 
tant but sometimes neglected effect when testing for the 
presence of a SIDM constant density core. 

Direct constraints from uncorrelated data. We conduct 
minimal data processing. We use the mass models to 
calculate, project, and PSF-distort theoretical spectra, 
and then fit these to the measured X-ray count spectra. 
We avoid deprojecting the data, thus guaranteeing fits to 
uncorrelated data. The gas temperature is handled as an 
internal variable, so that we skip the temperature fitting 
stage altogether. For Sunyaev-Zel'dovich measurements, 
we calculate and fit the uncorrelated Fourier modes of 
the decrement directly. 

No X-ray temperature weighting. Because the mass 
models arc fit directly to the spectra there is no need to 
calculate emission- or otherwise-weighted temperatures. 
Hydrodynamic N-body simulations show that emission- 
weighted temperatures do not accurately reflect the mea- 
sured spectroscopic temperature (Mazzotta et al. 2004; 
Rasia et al. 2005); recent methods for getting around 
this problem require the calculation of theoretical weights 
that depend on the X-ray calibration files (Vikhlinin 
2006). JACO calculates 2D spectra directly from the 
mass models. Temperature and density variations within 
any given annulus are reflected in the projected 2D X-ray 
spectra without any emission weighting. 

What is truly new here is (1) the direct calculation of 
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X-ray spectra from the gravitational potential, (2) the 
use of this potential to model Simyaev-Zel'dovich and 
weak lensing observations simultaneously, and (3) the 
ability to conduct a full, covariant error analysis on all 
the physical parameters. There are also other important 
features that have not been implemented in any prior 
method, e.g. the strict calculation of the mean molec- 
ular weight from the metallicity profile (§2.2) and the 
formal modeling for the termination shock as the clus- 
ter boundary (§2.3). Some of the other aspects of JACO 
have been discussed before. For example, Markevitch & 
Vikhlinin (1997) separate the gravitational potential into 
gaseous and dark components when modeling Abcll 2256. 
Fukazawa et al. (2006) and Humphrey ct al. (2006) em- 
ploy the stellar profile of the BCG in the X-ray fit. And 
"Smaug" (Pizzolato et al. 2003) was the first technique 
to incorporate direct fitting of projected spectra. JACO 
improves on Smaug by fitting a complete physical model, 
and avoiding the use of emission-weighted temperatures 
and coarsely gridded cooling functions. While JACO in- 
corporates relevant ideas from its predecessors, it was 
written completely independently, as the physical mo- 
tivation behind JACO made reuse of other code either 
slow or impossible. 

2.2. Fundamental Equations 
We first write the emissivity of the intracluster plasma, 

e v =n e n H \ u {T, Z) = ( — ) n 2 H \ u {T, Z), (1) 
\n H / 

where T is the temperature, n e is the electron density, 
and Uh is the hydrogen nucleus density. The cooling 
function A„ (T, Z) is weighted over all metals with rela- 
tive abundances fixed on the Grevesse & Sauval (1998) 
scale; in this system the absolute abundance is Z. For 
X v , JACO users can choose between the APEC and the 
MEKAL plasma codes, or they can provide their own. 
Note that in ionization equilibrium {n e /nu) is a function 
of metallicity; JACO recalculates it each time it derives 
an X-ray spectrum for a new value of Z. The same ap- 
plies to the mean molecular weight p. We find that if we 
did not recalculate p and (ne/nn) each time (i.e., held 
them fixed at some fiducial Z), we would be making up 
to a 5% systematic error in the output spectrum. 

In the approximation that the cluster contains a re- 
laxed, spherical, and stationary ideal gas, 
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where m p is the proton mass, Md is the dark mass, M g is 
the gas mass, and M s is the stellar mass contained with 
a radius r. Axisymmetric solutions will appear in future 
work. 

The most general solution to equation (2) is 
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where P c is the pressure at an arbitrary radius r = r c . 
The emitted spectrum within any volume then becomes 
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dV (4) 



Thus the spectrum is expressed purely as a function of 
mass and metallicity. No temperature fitting is required. 

Once we know the pressure, the Sunyaev-Z'eldovich 
decrement follows directly from an integral of the pres- 
sure along the line of sight (Birkinshaw 1999): 

Vsz = [ —^jdz (5) 

where <jt is the electron scattering cross section, and m e 
is the electron mass. Similarly, the reduced weak gravi- 
tational lensing shear in the tangential direction may be 
calculated directly from the total mass model (Miralda- 
Escude 1991): 
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AttG D A p WL 

Here qt is the reduced tangential shear, k is the con- 
vergence, S(i?) is the surface mass density, Da is the 
angular diameter distance to the cluster, and Pwl is a 
measure of the redshift distribution of the lensed sources; 
see Hoekstra et al. (1998) (§6), and §3.3 below. 

2.3. Boundary conditions 

The values of P c and r c arc influenced by the choice 
of the outer cluster boundary. In the idealized case that 
the cluster is infinite, P c and r c are fixed by the physical 
requirements that the temperature be everywhere finite 
and that the gas density decline monotonically with r. 
Then as r — > oo, the term in parenthesis must vanish at 
least as fast as p g ; otherwise, the temperature at large 
radii will be infinite. Thus 

GMp Ur> (10) 



and 



kT(r) = 



jj,m p 



GMpg 



dr'. 



(11) 



Of course, real relaxed clusters are not infinite; at 
some point, accretion from the surrounding intergalac- 
tic medium, and its associated shocks, become impor- 
tant (Ostriker, Bode, & Babul 2005). To simulate this 
we truncate the cluster at its virial radius. P c is then 
the "surface pressure" or pressure of the gas at the outer 
boundary. Assuming that the density of the intergalactic 
medium and the ICM at the termination point are the 
same, a useful estimate of this boundary pressure is 
p 2 GM tot 

Pc = qPgV cu 



QPg- 



:i2) 



Where r v i r is the virial radius (and the radius of termina- 
tion), i; c i rc is the circular velocity of the halo, and q is a 
constant of order unity. Both q and function of 

the cosmological parameters. For the prevailing ACDM 
cosmology, q ~ 1.1 and r v j r = rioo 5 the radius at which 
the cluster density is 100 times the critical density of the 
universe (Pierpaoli, Scott, & White 2001). 
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2.4. Model Radial Profiles 

We now require parameterized profiles to insert into 
the right-hand sides of equations (1) and (3). Because we 
are eliminating the temperature from the two equations, 
no assumptions regarding the form of the temperature 
profile arc required. However, we still require parame- 
terized gas, dark, and stellar mass profiles, as well as a 
metallicity profile. 

For the gas profile, we use a "triple" /3-model gas den- 
sity profile. Most relaxed clusters of galaxies studied 
to date have surface brightness distributions that are 
well-described by a double /3-model (Mohr, Mathiesen, 
& Evrard 1999; Hicks et al. 2002; Lewis et al. 2003; 
Jia et al. 2004; Johnstone et al. 2005); however certain 
cases a triple version is required to achieve a good fit 
(Pointccoutcau ct al. 2004; Vikhlinin et al. 2006, e.g.). 

Our version of the triple /3-model is 



■3ft /2 



(13) 



Note that in contrast to previous work, we define the 
triple (3 model in density, not in surface brightness. The 
asymptotic behavior of our triple (3 model is the same as 
the behavior of a surface brightness model with the same 
(3 parameters; only the details of the transitions among 
the various regimes are different. While the multiple (3 
model has several more fitting parameters than a broken 
power law, it is demonstrably better in many cases, e.g. 
Sun et al. (2004). 

For the dark matter distribution, JACO allows for a 
choice of profiles from the N-body and dynamical lit- 
erature. A fundamental prediction of most collisionless 
CDM simulations is that pdmM rises inwardly as or 
steeper to radii comparable within the resolution limit. 
The resulting model has become known as the "univer- 
sal" dark matter profile: 



Pu = Po{r/r d ) ™(l + r/r d ) 



(14) 



where x = r/r d . Traditionally, n = 1 (Navarro et al. 
1997) to 1.5 (Moore et al. 1998); however, we explore 
all values < n < 2 to test for softened as well as cuspy 
dark matter potentials. A similar profile that closely 
describes the distribution of stars in elliptical galaxies is 
the "7 profile" (Dehnen 1993; Tremaine et al. 1994): 



P 7 =Po(r/r d )-"(l + r/r d )"- 



(15) 



The widely used Hcrnquist (1990) profile is a special case 
of the 7 profile with n = 1. Unlike the universal profiles, 
the 7 family of profiles have simple analytical properties 
and finite total mass. We also allow for a simple, single 
power law mass distribution by taking the limit r d — > oo 
in equation 14. 

Finally, the stellar matter distribution in the BCG can 
be modeled either as a single /3-model or as a modified 
three-dimensional Sersic (1968) profile: 



p s = p cxp-(r/r s ) c 



(16) 



Thus a s = l/n s , where n s is the Sersic index such that 
n = 4 gives a de Vaucoulcurs (1953) profile. The shape 
parameters of the this profile are measured from high 
quality optical photometry, while the normalization is 
allowed to vary freely in the fit. The BCG profile allows 



us to define a constant baryon fraction (/&) dark matter 
model — this is a fiducial model where the dark matter 
density is proportional to the sum of the gas and the 
BCG stellar densities. 

The 7 and Sersic profiles above have closed-form inte- 
grated mass profiles. For the others, the enclosed mass 
can be expressed as a rapidly computable incomplete 
beta function (see Appendix A). 

The radial dependence of the metallicity distribution is 
not as well constrained by previous studies. We do know 
that relaxed clusters of galaxies show ample evidence of 
negative radial metallicity gradients, but the exact form 
of the decline is not clear. We choose the general distri- 
bution 

, , *<" = <"> 

i.e., the metallicity progresses smoothly from Zq at the 
center to Z\ at the outer regions, with a rate controlled 
by rz- This functional form is consistent with the abun- 
dance profiles derived from observations (Irwin & Breg- 
man 2001; De Grandi & Molendi 2001; Sun et al. 2003; 
Buote et al. 2003; de Plaa et al. 2004). 

Each radial profile has 3 or more free parameters. The 
grand total fit to the cluster involves more than a dozen 
free parameters. Table 1 is an annotated list of all such 
parameters. Not all parameters need to be fit at the same 
time; any number can be linked together or frozen at 
specific values. In most applications, the cluster redshift 
and the parameters relating to the shape of the stellar 
light profile will be fixed at specific values, because it is 
difficult or impossible to estimate them from X-ray data 
alone. In the case of poor quality data one can restrict 
the cluster model to a single /3-model or to a constant 
metallicity. 

2.5. Projection 

Dcprojected spectral or temperature data points are 
always correlated. Use of the x 2 statistic to fit correlated 
data is complex, and problematic if the covariance matrix 
of the deprojected data is unknown. A straightforward 
alternative is to project the model spectra on the sky, 
and fit them to the uncorrelated data. For projection we 
follow the "Smaug" method (Pizzolato et al. 2003) in 
detail. If the inner and outer edges of the ith annulus 
are at projected distance Ri-i and Ri respectively, then 
the observed flux is: 

1 f Ri /''"max Of , r A r 

F * ] =^iL^ mR L 7w=w (18) 

Where v' = v/(l + z), Dl is the luminosity distance 
to a source at redshift z, and r max is the outer limit of 
the cluster. By switching the order of integration, it is 
possible to further simplify the integral (Pizzolato et al. 
2003): 1 



Mi) 



D\ ! Ri _ 



r e v iK(r)dr 



where 



K(r) 



r 2 -RU 



ifr < Ri 



yV 2 - - v / r 2 - Kj ifr > R 



(19) 



(20) 



1 Note that Pizzolato et al. (2003) have an extra 47T in front of 
their equation (3) because they take a different definition of the 
cmissivity. 
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TABLE 1 

List of Parameters in the Cluster Model 



Parameter 


Sherpa Name 


Units 


Can be fit? 


Comments 


A 


contrast 




N 


Overdensity for masses and concentrations 




precision 




N 


Precision goal 


Pd 


darkmodcl 




N 


Integer specifying dark mass model 


p* 


starmodel 




N 


Integer specifying stellar mass model 


z 


redshift 




Y 


Cluster redshift 


Mg 


Mg 


10 14 M Q 


Y 


Gas mass within 




rxl 


Mpc 


Y 


Core radius of first /3-model 




rx2 


Mpc 


Y 


Core radius of second /3-modcl 


r X3 


rx3 


Mpc 


Y 


Core radius of third /3-model 


Pi 


bl 




Y 


Slope of first /3-modcl 


fa 


b2 




Y 


Slope of second /3-model 


fa 


b3 




Y 


Slope of third /3-modcl 


h 


xfrac2 




Y 


Normalization of second /3-model relative to the 1st 


f 

73 


xfrac3 




Y 


Normalization of third /3-model relative to the 1st 


M d 


Md 


10 14 M Q 


Y 


Dark mass within 


n 


ndark 


Y 


Slope of dark matter density profile 


ro 


rO 


Mpc 


Y 


Size of inner constant-density dark matter core 
Dark matter concentration 21 




rdml 




Y 


Z 


zO 




Y 


Mctallicity at r = 


Zoo 


zinf 


10 14 M Q 


Y 


Mctallicity at large radii 


Mbcg 


Mstar 


Y 


Total stellar mass of BCG 




alpha 


Y 


Index of 3D stellar light profile, cxp —(r/rs) as 


r s 


rs 


Mpc 


Y 


Scale radius of stellar profile 


^"shock 


rshock 


Mpc 


Y 


Radius of termination shock 


Note. - 


Additional parameters 


representing the residual soft X-ray background may also be fit to the data. 


a The user 


can decide to fit the concentration c = r^/r^, 


or else to fit r d itself. 





2.6. PSF Distortion 

A crucial step in comparing the projected model to 
the data is to account for PSF distortions. The JACO 
approach is to manipulate the data as little as possible, 
so we apply the distortion to the model spectra F v (i) via 
a convolution matrix: 

N 

S v (i) = J2^(iJ)FAj) (21) 
j=i 

where S v (i) is the final model spectrum for compar- 
ison with the data, and U u (i,j) contains the energy- 
dependent contribution of each annulus to itself and ev- 
ery other annulus. Given any set of annuli it is possible 
to calculate H u (i,j) and thus convolve the model with 
the PSF. 

In general, the PSF of Chandra and XMM-Newton is 
a function of energy and position. JACO makes an al- 
lowance for this. However, the shape of the PSF is not 
strictly circular or even ellipsoidal. We take no account of 
this anisotropy, instead treating the PSF as azimuthally 
symmetric at every point. In-flight XMM-Newton cal- 
ibration studies have found that this treatment of the 
PSF yields sufficiently accurate encircled energy profiles 
except for very bright point sources (Ghizzardi 2001a, 
2001b). The model PSF is given by 

p„(0) cx (1 + 2 /0o) C (22) 

where <j> is the angular distance in arcseconds between 
the source center and the position being considered. The 
shape parameters (fio and £ are functions of energy and 
position across the detector, and are empirically fit as 
polynomials in photon energy E = hv and source posi- 
tion 9: 

(pa = ci + c 2 E + c 3 9 + c A E9 (23) 



TABLE 2 
Energy Dependence of the PSF 



Coefficient 


ACIS 


MOS1 


MOS2 


pn 


ci 


1.137 


5.074 


4.759 


6.636 


C2 


-0.054 


-0.236 


-0.203 


-0.305 


C3 


0.076 


0.002 


0.014 


-0.175 


C4 


0.034 


-0.018 


-0.023 


-0.007 


C5 


6.119 


1.472 


1.411 


1.525 


C6 


-0.254 


-0.010 


-0.005 


-0.015 


C7 


-0.652 


-0.001 


-0.001 


-0.012 


C8 


0.051 


-0.002 


0.000 


-0.001 



C = c 5 + c 6 E + c 7 6 + c 8 E0 (24) 

Here 9 is measured in arcminutes and E is measured 
in keV. The XMM-Newton calibration team has derived 
(Ghizzardi 2001a, 2001b) values of ci_ 8 from archival ob- 
servational data. While extensive treatment of the Chan- 
dra PSF exists, no functional fits to the azimuthally av- 
eraged Chandra PSF have been published. For this rea- 
son, we undertake ray-tracing MARX (Wise 1997) sim- 
ulations of bright point sources at various off-axis angles 
and energies. 

We find that equation (22) provides an acceptable de- 
scription of the Chandra PSF as well. We measure ci_§ 
in a similar manner to Ghizzardi (2001a), except that we 
use ray-traced rather than in-orbit data. The polynomial 
coefficients for all instruments are shown in Table 2. 

Using this functional form for the PSF we can calculate 
the PSF convolution matrix Tl v (i,j). The details of the 
calculation are given in Appendix B. 

2.7. Software Interface 
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0.01 0.1 

Projected Distance (arcmin) 

Fig. 1. — (Zeft) HST WFPC2 image of the core of Abell 478, with superposed X-ray contours from the Chandra ACIS-S image. 
Profile of the surface brightness plus background for the BCG, as well as the best-fit Sersic model (a = 0.289). 

The JACO software package will soon be publicly avail- 
able at http://jacocluster.sourceforge.net. It consists of 
the following components: 

1. A core C language library which calculates the 
observed X-ray spectrum, tangential shear, and 
Compton y parameter from the input dark, gas, 
and stellar mass distributions; 

2. An interface which links the core library to the 
Shcrpa data analysis package (Freeman, Doe, & 
Siemiginowska 2001) as a standard user-defined 
model, and provides the necessary scripts for 
graphical viewing and fitting; 

3. A parallelized interface, Hrothgar, which runs 
the core JACO routine on multiprocessor facili- 
ties such as Beowulf clusters. Hrothgar is a gen- 
eral purpose package and will soon be available at 
http : / /hrothgar . sourceforge .net . 

4. A set of optional data reduction scripts which pro- 
cess standard archival X-ray data releases accord- 
ing to the procedures described in Appendix C. 
Weak lensing and SZ data reduction is up to the 
user. On request, the authors can also provide a 
routine for transforming JACO Compton y maps 
into interferometer observables. The existing rou- 
tines are customized for the Cosmic Background 
Imager (Padin et al. 2002), but other arrays are 
easily accommodated. 

JACO offers a choice of numerical integration meth- 
ods. The faster method uses an adaptive Gauss-Legendre 
quadrature with 10 or 20 abscissae depending on the 
complexity of the integral. JACO also offers adaptive 
Gauss-Kronod quadrature with convergence checking. 
The latter method is slower but is guaranteed to con- 
verge to a given accuracy. We find that the faster method 
produces model spectra with a median accuracy of 0.1%. 

3. APPLICATION TO ABELL 478 

Because the JACO technique is new, we apply it to 
a well-known, relaxed cluster of galaxies. Abell 478 is 



right 



a rich cluster at z=0.088. The most recent available 
studies with Chandra (Sanderson et al. 2005; Voigt & 
Fabian 2006) and XMM-Newton (Pointecouteau et al. 
2004) suggest a peak ICM temperature of of ~ 7 — 9 
keV, making this cluster one of the more massive known 
within z < 0.1. For this analysis we assume Hq = 70 km 
s^ 1 Mpc -1 , tt = 0.3, and Q A = 0.7. With these values, 
1' = 98.7 kpc at the distance of Abell 478. 

This cluster is an ideal test case for JACO, and high- 
lights the advantages of direct spectral fitting for hydro- 
static mass determination. The Chandra and XMM- 
Newton observations record a total of over 1,000,000 
source photons; yet Pointecouteau et al. (2004), Sander- 
son et al. (2005), and Voigt & Fabian (2006) show that 
temperature profiles of this cluster disagree when calcu- 
lated with single temperature emission-weighted plasma. 
We will show that the addition of optical data (for the 
stellar mass profile and the weak gravitational lensing 
shear) and radio data (for the Sunyaev-Zel'dovich decre- 
ment) improves the constraints on the dark matter distri- 
bution in the cluster, and bring the Chandra and XMM- 
Newton data into closer agreement. 

3.1. X-ray Data 

To model Abell 478, we fit JACO X-ray models to 
available XMM-Newton and Chandra archival. There 
are four instruments in all: EPIC pn, EPIC MOS 1 and 
2, and Chandra ACIS-S. We used XMM-Newton obser- 
vation sequence 109880101 (126ks, of which 43ks were 
useful) and Chandra ObsID 1669 (42ks, of which all was 
useful). For the Chandra data, we use CIAO 3.3 and 
CALDB 3.2.1. For the XMM-Newton data, we use the 
latest CCF calibration release before August 2006. 

We apply the generic reduction described in Appendix 
C. We extract spectra in circular regions around the 
X-ray centroid, which coincides in all datasets with 
RA(J2000) = 04:13:25.3, DEC(J2000) = +10:27:54. A 
total of 76 annuli are used, of which Chandra covers only 
the innermost 40. The ACIS data are used to choose the 
sizes of the innermost 40 annuli; we require at least 2500 
background-subtracted counts per annulus. A similar 
requirement is used to choose the sizes of the outer 36 
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Fig. 2. — (left) SZ data and residuals for the best-fit JACO model; (right) weak lensing data and residuals for the best-fit JACO model 



annuli -2500 background-subtracted counts in the EPIC 
pn camera. The innermost annulus is 13" (22 kpc) in 
radius; the annuli increase regular by 5" (8 kpc) until a 
distance of 4.6'(450 kpc), at which point they increase by 
14"-40" until the final annulus at 9' (890 kpc). The ex- 
tracted spectra were fit over the 0.6-6 keV energy range. 
The JACO plasma code was set to MEKAL. 

The galactic absorption column varies with position 
on the sky (Sun et al. 2003; Pointecouteau et al. 2004; 
Vikhlinin et al. 2005). This effect is more important for 
the XMM-Newton data, which cover a wider field than 
Chandra; we take the gradient into account by allowing 
the hydrogen column density to vary linearly with radius. 
Taking the variation of nn with position into account still 
does not resolve the temperature discrepancies between 
the Chandra and XMM-Newton temperatures at inter- 
mediate radii (Pointecouteau ct al. 2004; Vikhlinin et al. 
2005). 

3.2. Hubble Space Telescope Data 

Figure 1 shows a 300s archival Hubble Space Tele- 
scope WPFC2 image of the central <~ 60 kpc diameter 
region. The corresponding HST observation number is 
U5A40901R, taken with the - 20001 bandpass F606W 
filter. Superposed are the Chandra ACIS surface bright- 
ness contours. 

The stellar mass of the BCG makes a non-negligible 
contribution to the matter distribution within 50 kpc 
(Sand et al. 2002, 2004; Vikhlinin et al. 2006), and hence 
contributes to the equation of hydrostatic equilibrium at 
all radii. To model the stellar mass profile of the BCG, 
we fit its HST surface brightness profile with a 3D Sersic 
(1968) model (equation 16). Figure 1 shows the best 
fit model: a BCG = 0.289 ± 0.002 and r = 0.36 ± 0.02". 
This model is added to the gaseous and dark components 
when conducting the grand total fit. 

Using the photometric solution for the HST image, we 
calculate an uncorrected F606W magnitude of 15.1 ±0.1 



by integrating the best-fit Sersic model to infinity. The 
galaxy is in a high-extinction region; E(B-V) = 0.589 
(Schlegcl, Finkbeiner, & Davis 1998). After applying 
the extinction correction of 1.7 mag and a k-correction 
0.14 mag (Poggianti 1997), we find that the total F606W 
luminosity of the galaxy is 4.9 ± 0.3 x 10 11 Lq. This 
value is important in constraining the average stellar 
mass-to- light ratio of the BCG, Tg CG . For an evolved 
stellar population with mean age > 3Gyr in the wave- 
length regime of the F606W filter, it is expected that 
1 < Tg CG < 4Mq/L depending on the stellar initial 
mass function (IMF) (Maraston 1998). The lower limit 
of 1 applies regardless of the choice of IMF (Maraston 
1998). 

3.3. Weak Lensing Data 

The weak lensing measurements for A478 are based 
on archival data taken with the CFH12k camera on the 
Canada- France-Hawaii- Telescope (CFHT). The camera 
consists of an array of 6 by 2 CCDs, each 2048 by 4096 
pixels. The pixel scale is 0.206'.', which ensures good 
sampling for the sub-arcsecond imaging data used here. 
The resulting field of view is about 42 by 28 arcminutes, 
which is of great importance when studying nearby clus- 
ters, such as A478. 

The weak lensing analysis requires good image quality 
and therefore we selected only data with seeing better 
than 0.8". This selection yielded 9 exposures with 605s 
of integration each, resulting in a total integration time 
of 5445s. Unfortunately only i?-band observations were 
available. This in principle complicates removal of clus- 
ter members. However, our study of a large suite of R- 
band observations of other clusters (Hockstra 2007) has 
shown that we can readily correct for this source of con- 
tamination using the excess galaxy counts as a function 
of radius. 

Detrended data (de-biased and flatfielded) are pro- 
vided to the community through CADC. We process this 
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Fig. 3. — Spectra and 
the 268 total spectra are 



best-fit models for Chandra ACIS-S (as), XMM-Newton MOS1 (ml), MOS2 (m2), and pn cameras. Only 72 of 
shown; the inner and outer radii of the extracted annuli (in arcminutcs) appear next to the instrument name. 



data through the analysis pipeline described in Hoekstra 
et al. (1998), Hoekstra, Franx, & Kuijken (2000), and 
Hoekstra (2007). First we use the hierarchical peak find- 
ing algorithm from Kaiser, Squires, & Broadhurst (1995) 
to find objects with a significance > 5a over the local sky. 
The peak finder gives estimates for the object size, and 
we reject all objects smaller than the size of the PSF. The 
remaining objects arc analyzed, which yields estimates 



for the size, apparent magnitude and shape parameters 
and their measurement errors. The image is inspected by 
eye, in order to remove areas that would lead to spurious 
detections. 

To measure the small, lensing induced distortions it 
is important to accurately correct the shapes for PSF 
anisotropy, as well as for the diluting effect of seeing. To 
characterize the spatial variation of the PSF we select a 



JACO 



9 




R (Mpc) 

Fig. 4. — Two visualizations of the best-fit universal models (equation 14) for Abell 478. (top) The spectra are collapsed into three color 
regimes: 0.6-1.5 keV (upper line), 1.5-2.0 keV (middle line), and 2.0-6.0 keV (lower line). The count rates are offset 0, -1, and -2 orders of 
magnitude, respectively, for clarity. Data to model ratios are also shown, (bottom) The spectra and models are shown as color ratios. We 
show the ratio of count rates in the three bands to the total count rate: 0.6-1.5 keV (lower line), 1.5-2.0 keV (middle line), and 2.0-6.0 keV 
(upper line). Each error bar here represents ~ 50-70 of the fit spectral data points, each of which in turn represents ~ 50 photons. 



sample of moderately bright stars from our images (e.g. 
Hoekstra et al. 1998). Seeing circularizes the images, 
thus lowering the raw lensing signal. To correct for the 
seeing, we need to rescale the polarizations to their "pre- 
seeing" value, as outlined in Hoekstra et al. (1998). Our 
pipeline has been tested on simulated images in great 
detail (e.g. Hoekstra et al. 1998; Heymans et al. 2006). 
These results suggest that we can recover the shear with 



an accuracy of ~ 2% (Heymans et al. 2006). 

The catalog of objects with corrected shapes is used 
for the weak lensing analysis. In the analysis presented 
here, we consider the tangential distortion as a func- 
tion of radius from the cluster center. The resulting 
measurements, using galaxies with apparent magnitudes 
22 < R < 24.5 are shown in Figure 2. 

The interpretation of the lensing signal requires knowl- 
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edge of the rcdshifts of the source galaxies. Based on our 
data alone, we do not know the redshifts of the individ- 
ual sources. The observed lensing signal, however, is an 
ensemble average of many different galaxies, each with 
their own redshift. As a result, it is sufficient to know 
the source redshift distribution to compute the average 
value (3 = (Di s /D s ). Based on the Hubble Deep Fields, 
we obtain (3 = 0.8. We note that, because A478 is a 
low redshift cluster, the inferred lensing mass is rather 
insensitive to the detailed source redshift distribution as 
effectively all background galaxies are far away. 

3.4. Sunyaev-Zel'dovich Data 

The SZ data were taken with the Cosmic Background 
Imager (CBI; Padin et al. 2002) over 11 nights as part of 
a complete sample of X-ray luminous, low-redshift (z < 
0.1) clusters. A478 is a remarkably clean cluster, with 
no apparent point sources at 30 GHz in the CBI data. 
By far the dominant source of noise in the data is the 
CMB, which reduces the overall significance of detection 
from 24.8 sigma, if only thermal noise is considered, to 
8.3 sigma. See Udomprasert et al. (2004) for a more 
complete treatment of observing the SZ effect with the 
CBI, as well as a more detailed description of the A478 
data used here. 

To calculate x 2 from the SZ data, we first run the 
visibilities through CBIGRIDR, the CBI CMB analy- 
sis pipeline (Myers et al. 2003). This compresses them 
into 2000 " gridded estimators" , and calculates the noise 
and CMB covariances between the estimators, as well 
as source vectors for source projection if so desired (e.g. 
Bond, Jaffc, & Knox 1998). This is done once per cluster, 
before fitting cluster models. During the fit, Compton 
y maps for models are run through the CBI simulation 
pipeline to generate predicted visibilities, which arc then 
run through CBIGRIDR to generate the predicted esti- 
mators. Let A be the estimators, Am be the predicted 
estimators for the current model, Scm b be the CMB co- 
variance matrix, and N be the noise covariance matrix, 
then x 2 = (A - A M ) T (N + Scmb)~ 1 {A - A M ). Since 
the noise and the CMB covariance are independent of 
the cluster model, the inverse need be taken only once. 
In practice, we rotate into a space in which N + Scmb is 
diagonal using V, the eigenvectors of N + Scmb- Let 
A* = V T A, then X 2 reduces to £(A* - A Ml ) 2 /af 
where of is the corresponding eigenvalue of A + Scmb- 
This lets us treat the rotated estimators as independent, 
uncorrelated measurements of the sky. In practice, it 
takes a few minutes to calculate the eigenvector rotation, 
and each x 2 evaluation (including running the visibility 
pipeline and CBIGRIDR for estimators) takes ~ 0.5 sec- 
onds. The data for the mode numbers which contain 
most of the cluster signal appear in Figure 2. We do not 
fit mode numbers < 600, which are guaranteed not to 
contain any cluster signal. 

3.5. Results 

First we fit the JACO models separately to the Chan- 
dra and XMM-Newton data. We thereby gain insight 
into how the differences among the four separate instru- 
ments affect the estimated physical parameters. The 
breakdown for each instrument as well as the number 
of radial bins appear in Table 3. We show ~ 25% of 



the X-ray spectra along with the best-fit spectral mod- 
els in Figure 3. We also show collapsed views of all 
the spectral profiles in Figure 4. The number of radial 
bins exceeds those in previous mass measurement pa- 
pers (Pointecouteau et al. 2004; Voigt & Fabian 2006; 
Vikhlinin et al. 2006) by a factor of <~ 5. This allows 
us to constrain the dark mass and gas mass (or surface 
brightness and temperature) profiles simultaneously, and 
ensures that we do not miss any relevant details in the 
spatial distribution of the spectra. The simultaneous fits 
to the SZ and weak lensing data are shown in Figure 2. 

The constraints from the fits appear in Figure 5, where 
we show 68% and 95% the confidence intervals for some 
of the chief parameters of interest. All plots assume a 
universal profile with free inner slope n (equation 14). 
We show the constraints both with and without the in- 
clusion of the weak lensing and SZ data. 

The differences in the calibration among the various 
four instruments are immediately apparent in the X-ray 
fits. For example, as is typical of rich, hot clusters, the 
Chandra data yield a higher average temperature (Ko- 
tov & Vikhlinin 2005; Vikhlinin et al. 2005), resulting 
in a higher dark mass measurement. Also notable is the 
fact that the errors in the Chandra-derived quantities are 
larger. Even though the Chandra ACIS-S has the supe- 
rior spatial resolution, XMM-Newton has the greater sky 
area, so that the constraints on the shape of the dark pro- 
file, and especially on the masses at f2500j are tighter for 
the XMM-Newton data. This is due to the fact that our 
outer extraction radius for Chandra (450 kpc) is at about 
0-6r 2 5oo; quantities computed outside this radius are ex- 
trapolations and therefore subject to greater uncertainty. 

It is also clear from Figure 5 that the addition of the 
lensing and SZ data can help constrain the dark matter 
parameters significantly. For example, the uncertainty 
in the dark matter concentration as measured by Chan- 
dra is halved through the addition of the lensing and SZ 
data. We explore the additional power afforded by these 
types of observations in §4.4 below. Most encouragingly, 
there is no bias or disagreement whatsoever in the joint 
fit with the additional data sources — the same physical 
model can account for the X-ray, SZ, and weak lensing 
observations at the same time. 

Chandra is the instrument for which n is most affected 
by the addition of the SZ and lensing data. The dark 
matter slope as measured by Chandra alone (n < 1) 
is consistent with the XMM-Newton value (1.1 ± 0.3). 
When we add the SZ and lensing data to the Chandra 
data, the slope (n = 0.7 ± 0.4) is brought into closer 
agreement with XMM-Newton, in that the low n = 
solution is disfavored. 

4. DISCUSSION 

4.1. Comparison with Previous X-ray Studies 

For the first time, we are able to simultaneously con- 
strain all the parameters of physical interest in the mass 
fitting process, some of which are listed in Table 3. JACO 
allows us to calculate the joint probabilities of quantities 
such as the dark mass, the gas metallicity, and the ab- 
sorbing hydrogen column density. Determination of the 
covariance of such parameters is a unique property of the 
method, and cannot be easily be duplicated with previ- 
ous techniques. 

In Table 4 we compare the chief structural properties 
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Fig. 5. — 68% and 95% confidence contours on the dark mass within r2500i the gas mass within T2500, the smallest of the three /3 
parameters f}\, the concentration with respect to rsoo, the slope of the dark potential n, and the stellar mass-to-light ratio of the BCG 
T BCG- AU masses are in unit s of 1O 14 M . Clockwise from the top left are the Chandra ACIS-S, XMM-Newton MOS1, XMM-Newton pn, 
and XMM-Newton MOS2 instruments. The filled contours show joint X-ray, SZ, and weak lensing constraints. The unfilled contours show 
the X-ray constraints alone. The thick lines show the allowed values (from theory) for Tg CG . 



of the best-fit JACO model with previous studies of Abell 
478. In general, the JACO analysis of the Chandra and 
XMM-Newton data is in excellent agreement with pre- 
vious works. In comparing the JACO Chandra results 
with Sanderson et al. (2005) and Vikhlinin et al. (2006), 
we note that the uncertainties in the measured JACO 
parameters are larger. The explanation for this is that 
JACO allows all physical quantities to vary simultane- 
ously during the fitting process. Other fitting techniques 
conduct a multi-stage fit: first they fix the surface bright- 
ness profile at the best-fit multiple (3 model, and then, 
using this fixed profile, they fit a temperature profile to 
the reduced spectra. This results in an underestimate 



of the true uncertainties in n and the mass, which are 
better reflected in the JACO measurements. 

We note that the JACO results isolate the dark mat- 
ter from the gas and stellar profiles. In other studies, 
the measured n reflects the total gravitational potential. 
The similarity of the JACO n values compared with the 
globally measured n in the literature yields the useful 
conclusion that the stellar and gas profiles do not signif- 
icantly affect the measured dark matter slopes, at least 
for Abell 478. 

4.2. Contamination by the Central Source 
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TABLE 3 

Best Fit Universal Models for Abell 478 



Chandra EPIC MOS1 EPIC MOS2 EPIC pn Joint Constraints 

X-ray x 2 /DOF 5294/5413 (g=0.87) 4170/4084 (<?=0.17) 6582/6326 (g=0.012) 7576/7462 (g=0.17) 23264/23345 (g=0.63) 

Joint x 2 / D OF 5659/5729 (g=0.74) 4526/4400 (<?=0.090) 6936/6642 (g=0.006) 7929/7778 (9=0.11) 23619/23661 (g=0.57) 

M g 25oo 0.403 ± 0.020 0.434 ± 0.013 0.414 ± 0.011 0.399 ± 0.009 0.425 ± 0.007 

M d 2500 4.273 ± 0.753 3.254 ± 0.263 2.748 ± 0.199 2.839 ± 0.164 3.181 ± 0.135 

n ' 0.709 ± 0.436 1.035 ± 0.312 0.665 ± 0.346 1.096 ± 0.322 1.012 ± 0.207 

C25oo 2.027 ± 1.270 0.836 ± 0.528 1.966 ± 0.893 1.045 ± 0.633 1.234 ± 0.489 

Z 0.931 ± 0.153 1.960 ± 0.474 0.731 ± 0.228 1.440 ± 0.718 1.082 ± 0.416 

Zoo 0.095 ± 0.488 0.350 ± 0.144 0.159 ± 0.114 0.443 ± 0.069 0.336 ± 0.065 

Tbcg < 3 2.573 ± 2.217 5.487 ± 1.829 4.280 ± 3.251 1.669 ± 1.390 



Note. — Joint constraints on parameters for the Universal fit to each datasct and to the joint data sets, including SZ and weak lensing 
observations. We also show the goodness of fit q for the X-ray observations alone. For a description of the parameters and their units, see 
Table 1 and §2.4. The confidence intervals were obtained via the fit covariance matrix. 



TABLE 4 

Comparison with Previous Work (X-ray Data Only) 





JACO 


S05 


VF06 


JACO 


P05 


JACO 


V06 


VF06 


JACO 


Parameters 3 - 


Chandra 


Chandra 


Chandra 


XMM 


XMM 


Chandra 


Chandra 


Chandra 


XMM 




All Data 


All Data 


All Data 


All Data 


All Data 


Excised 


Excised 


Excised 


Excised 


n h 


< 1 


0.35 ±0.22 


0.49 ± 0.45 


1.1 ±0.3 


= 1 


< 1.33 




1 1+ - 2 
i - i -0.6 


1.7±0.3 


^2500 


0.68 ±0.07 


0.62 ±0.07 


u - ' D -0.19 


0.59 ±0.02 




0.67 ±0.07 






0.60 ±0.03 


C2500 


3.4 ±2.1 




1.1 ±0.7 




< 10 






< 0.70 


M to t,2500 


4.9 ±0.9 






3.3 ±0.2 




4.7 ±0.8 


4.2 ±0.3 




3.4 ±0.3 


»*500 


1.5 ±0.2 






1.4 ±0.1 




1.4 ±0.3 


1.4 ±0.1 




1.5 ±0.5 


C500 


7.4 ±4.3 






2.4 ± 1.3 




8.6 ± 7.0 


3.8 ± 0.3 




< 2.0 


A^tot,500 


9.9 ±2.6 






8.2 ± 1.0 


7.6 ± 1.1 


9.2 ± 3.8 


7.8 ±1.4 




11±6 


»*200 


2.2 ±0.4 


2.2 ±0.1 


3.0+i 4 o 


2.1 ±0.2 


2.1 ±0.1 


2.1 ± 0.3 




Ri 9.2 


2.6 ±0.6 


C200 


11 ±6 


7±2 


2.91?/° 


3.6 ± 1.7 


4.2 ±0.4 


13 ±9 






< 1.2 


M to t,200 


13 ±4 


13 ±3 




12 ±2 


11 ±2 


12 ±3 






21 ±9 



Note. — Shown are 68% confidence intervals on the structural parameters. JACO: This work; S05: Sanderson et al. (2005); V06: 
Vikhlinin et al. (2006); VF06: Voigt & Fabian (2006); P05: Pointecouteau et al. (2005). "Chandra" or "XMM" refers to the observatory 
used by each of the cited works. "All Data" means that no portion of the cluster emission was removed during analysis; "Excised" means 
that the central as 20 — 30 kpc was removed before analysis. 

a Distances are in units of Mpc and masses are in units of 1O 14 M0. In all sources, measurements at rsoo and r200 are based on extrapolation 
from data within those radii. "Excised" fits exclude the central 30kpc. 

b Measurcd for the dark mass profile only in our work, and for the total mass profile in the other papers. All papers assume a universal 
profile (equation 14) with either free or fixed slope n; VF06 use a power law of slope n in their excised fit. 



A key problem with the mass measurement within the 
inner 20-30 kpc of Abell 478 (and other similar cool core 
clusters) is AGN activity. The symmetric bubbles dis- 
cussed by Sanderson et al. (2005) are clearly visible in 
Figure 1 . All X-ray mass measurement techniques rely on 
the assumption of hydrostatic equilibrium. However, in 
most relaxed clusters there is clear, complex interaction 
between the central AGN/radio source and the interven- 
ing gas (e.g. Blanton et al. 2003; Clarke, Blanton, & 
Sarazin 2004; McNamara et al. 2005). 

In the discussion thus far, we have conducted a hydrod- 
static analysis including this central region. There is at 
this point no consensus on whether hydrostatic analysis 
of noncquilibrium gas yields correct results. One com- 
mon approach (Voigt & Fabian 2006; Vikhlinin et al. 
2006) is to excise the inner region from analysis. Note 
that excising the X-ray emission around the BCG does 
not mean setting the mass within that region equal to 
zero — it merely means that the shape of the mass distri- 
bution inside that region is not constrained by the X-ray 
model. For XMM-Newton, excising only the central re- 



gion is problematic: the large size of the PSF and the 
steepness of the surface brightness profile ensure that 
the central spectrum contributes as far out as 1' from 
the center (Markevitch 2002). 

To examine the effects of removing the emission from 
the disturbed region, we repeat our entire X-ray analy- 
sis, excising the central 20" (32 kpc) for Chandra, and 
the central 1' (100 kpc) for the XMM-Newton data. The 
results appear in Table 4. We find that removing the cen- 
tral portion causes us to measure a steeper value for n 
in both the XMM-Newton and the Chandra data. Out- 
side 100 kpc, the XMM-Newton dark matter profile is 
consistent with a single power law (the concentration is 
consistent with in all cases) . The Chandra dark profile 
becomes more fully consistent with an n = 1 Universal 
profile, as first pointed out by Voigt & Fabian (2006). 

Given the quality of the Abell 478 observations, these 
results show that strong constraints on the shape of the 
dark profile depend greatly on the data within « 20 — 30 
kpc of the center. Unfortunately, we do not yet know 
whether the hydrostatic model considered here gives cor- 
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Fig. 6. — 68% and 95% confidence contours for the combined JACO fit to all X-ray, SZ, and weak lensing data. Shown are the dark mass 
within r2soo ( m units of 1O 14 M0), the gas mass within r2500i the three j3 parameters /3i,2,3, the concentration with respect to r2500i the 
slope of the dark potential n, the metallicity at the core (Zq) and at large radii (Zaa), the stellar mass-to-light ratio of the BCG Tg CG , the 
central galactic absorption column tihq in units of 10 22 cm -2 , and its gradient with respect to distance from the center 9 in arcminutes. 
The thick lines show the allowed values (from theory) for the stellar mass-to-light ratio Tg CG . 



rect results when applied to regions where AGN heating 
is important. However, at least for Abell 478, the mass 
profiles measured for the full data set are consistent with 
those measured with the central 30 kpc excluded (the lat- 
ter have larger errors). Thus, treating the central region 
as hydrostatic does not appear to bias the final measured 
parameters. 

4.3. Combined Constraints from All Data Sources 

We now consider the question of conducting a grand 
total fit to all four instruments. This is a problematic 
issue, because of the known disagreement in the temper- 



ature profiles of this cluster. In combination, the spec- 
tra from the four X-ray instruments, together with the 
lensing and SZ observations, contain 23345 data points 
(including over a million X-ray photons). To account for 
the differences in the cross-calibration, the simultaneous 
fit to all four instruments includes a 4% systematic error 
added in quadrature to the statistical error. Without the 
addition of a systematic component no model can simul- 
taneously provide a good fit to both the Chandra and 
XMM-Newton data. The magnitude of the systematic 
error was chosen by comparing count ratios in published 
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Fig. 7. — Testing the power of the SZ and weak lcnsing data 
without gas temperature information. Shown are the dark mass 
within 1 Mpc in units of 10 14 A/q, the dark radius from equation 
14 in units of Mpc, as well as the inner slope n. The X-ray gas mass 
is taken as a fixed prior, but X-ray temperature (spectral) informa- 
tion is not used. Solid contours show SZ+weak lensing constraints; 
dark unfilled contours show the lensing constraints alone, and the 
light unfilled contours show the SZ constraints alone. 

cross-calibration studies (Stuhlinger et al. 2006). 

Figure 6 shows the combination of all four instruments 
with the X-ray, SZ, and weak lensing data for the Uni- 
versal profile. The total number of degrees of freedom 
are 23661. We find that the n ~ 1 universal profile pro- 
vides the best overall fit, ^with \ 2 — 23619. The n ~ 1 
gamma-profile (essentially a Hernquist 1990, profile) pro- 
vides the next best match, with A% 2 = 15 with respect to 
the universal profile. The constant baryon fraction model 
and the single power law models provide worse matches 
to the data, with A\ 2 — 37 and A\ 2 = 187, respectively. 
Using the likelihood ratio test, and noting that a single 
power law is obtained from the universal profile by set- 
ting rd = 0, we estimate that the single power law model 
is disfavored with false-alarm probability p < 10~ 6 . 

4.4. The Leverage of SZ and Weak Lensing Data 

We showed in §3.5 that the addition of SZ and weak 
lensing data results in improved constraints on the shape 
parameters of the gas and dark matter density (at least 
for the Chandra data taken alone). Here we examine the 
reasons for this improvement in further detail. As an 
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Fig. 8. — Total enclosed mass for the dark, gaseous 
BCG components of Abell 478, in units of 10 14 M, 
the 68% confidence intervals for fit to all data (Chandra + XMM 
Newton + SZ + weak lensing). The stellar mass constraint includes 
the assumption that the stellar mass-to-light ratio Tg CG > 1. 

exercise, we fix the gas mass distribution at the profile 
derived from the fit to the combined X-ray data. We al- 
low only the dark and stellar mass profiles to vary. Thus 
we take advantage of the fact that the gas mass pro- 
file is a well measured and well agreed-upon quantity for 
most clusters. For this exercise we ignore the spectral 
information, where the Chandra and XMM-Newton dis- 
agreement is more serious. 

The results appear in Figure 7. These figures show that 
the SZ and weak Lensing data are complementary, as 
suggested in previous theoretical studies (Zaroubi et al. 
2001; Dore et al. 2001). While neither can by itself 
strongly constrain the slope and concentration of the 
dark matter profile, together they offer nearly orthog- 
onal constraints on these parameters. 

Comparing Figure 7 to Figure 5, it is clear that when- 
ever the quantity of X-ray data overwhelms the SZ and 
weak lensing data, the statistical weight of the latter two 
is small; this is why the high quality EPIC pn and MOS2 
results are nearly unchanged by the addition of the SZ 
and weak lensing information. However, when the spatial 
resolution and the number of X-ray photons is small — as 
is surely the case for most higher redshift clusters — the 
SZ and weak lensing information contribute substantially 
to the constraints on the dark matter profile. 

4.5. Dark Matter and the Role of the BCG 

The stellar mass of the BCG is an important consid- 
eration for modeling the inner regions of galaxy clusters 
(Lewis et al. 2003; Mamon & Lokas 2005). In Figure 
6, the correlation between the dark matter slope n and 
the stellar mass-to-light ratio Tg CG is strong. The more 
massive the BCG, the less room there is for a central 
dark matter cusp. On the other hand, for a given BCG 
luminosity, there is a minimum value of Tg CG regardless 
of the stellar initial mass function (IMF). The star light 
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in effect yields an upper limit on the steepness of the 
dark matter profile. 

A galaxy dominated by a > 3 Gyr old stellar popula- 
tion ought to have Tg CG w 1 - 4 in the HST F606W 
band (Maraston 1998). These limit are shown in Figures 
5 and 6. Adopting these limits yields further constraints 
on the dark matter slope n (towards lower values) and 
the concentration (towards higher values). For exam- 
ple, in Figure 6, values of n > 1.1 would be ruled out 
by the allowed Tg CG = 1 limit at the 99% confidence 
level. The ruled-out region changes to n > 1.3 if when 
we use the excised data only. It is worth noting that 
while the Maraston (1998) population synthesis models 
have different stellar mass-to-light ratios depending on 
the exact mixture of stellar populations, only the upper 
limit (Tg CG = 4) is sensitive to the shape of the initial 
mass function; the lower limit, Tg CG = 1, is firm and 
independent of the IMF for a given age. 

Figure 8 shows the mass fraction in stars for Abell 478 
as a function of distance from the center. As a contribu- 
tor to the equation of hydrostatic equilibrium, the stars 
in the BCG are as important or more important than the 
gas within 80 kpc. The stars also contribute 30% of the 
total mass within 20 kpc of the center. 

5. CONCLUSION 

We have developed a new method for measuring the 
dark matter profile of a cluster directly from X-ray, lens- 
ing, and SZ data. 

• JACO works directly in the data plane, and gen- 
erates the observables (X-ray spectra, weak lens- 
ing shear profiles, or Sunyaev-Zel'dovich decrement 
maps) from candidate models. It therefore allows 
simultaneous constraints on a cluster's dark matter 
profile using multiwavelength data. It allows joint 
constraints on all physical parameters of interest. 

• JACO explicitly separates the dark, gas, and stellar 
mass profiles. This allows extraction of the shape 
of the dark profile independently of the rest of the 
cluster. 

• As long as the gas mass profile is well known, the 
SZ and weak lensing data together can provide or- 
thogonal constraints on the shape parameters of 
the dark profile. 

• We present new CFHT weak lensing measurements 
for the well-studied rich cluster Abell 478, and 
provide an improved reduction of existing CBI 
Sunyaev-Zel'dovich data. We analyze these data 
in conjunction with existing Chandra and XMM- 
Newton observations of the cluster. We find ex- 
cellent agreement among all data sets when they 
are fit simultaneously with JACO models. The 



weak lensing and SZ observations improve the con- 
straints on the shape of the dark matter distribu- 
tion from the Chandra data. 

• The Chandra and XMM-Newton data for Abell 478 
yield similar slopes for the inner dark matter pro- 
file: n < 1 and n = 1.1 ± 0.3 at the 68% confidence 
level, respectively. The Chandra constraints be- 
come more fully consistent with the n = 1 NFW 
profile when we excise the morphologically dis- 
turbed central 30 kpc. A similar result was noted 
by Voigt & Fabian (2006). 

• At intermediate and higher redshifts, where the 
quality of the X-ray data rapidly diminishes, the 
SZ and weak lensing data become increasingly im- 
portant for characterizing the properties of dark 
matter. In these regimes JACO will be a powerful 
tool for improving constraints on the shape of the 
cluster potential. 

• JACO as described in this paper will 
soon be available for public download at 
http://jacocluster.sourceforge.net. 

In future papers in this series, we will test JACO 
against gasdynamical N-body simulations, generalize the 
method to the axisymmetric case, and use JACO to an- 
alyze multiwavelength data for a large cluster sample. 
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APPENDIX 



A. FAST EVALUATION OF BROKEN POWER LAW INTEGRALS 
To calculate the integrated mass of the profiles considered in §2.4, we often need to calculate integrals of the type 

f-r/r t 



M(r) 



f r f T ' rt n-q 

/ A-Kr 2 p(r/r t )dr = iirrfpo / x 2 ~ n (l + x p ) *> dx (Al) 
Jo Jo 
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where r t is the characteristic radius of the density p, and n and q are the slopes at r = and r = oo, respectively. 
Rather than resorting to hypergeometric functions to evaluate the integral, we can express it in terms of the incomplete 
beta function B y (a, b) for which very fast numerical routines exist: 

M(r) = 4nr* Po B y . ys (A2) 

In the common definition of the incomplete beta function, a and b must be positive definite. However, q — 3 could well 
be negative (e.g., for a /3-model). The following recurrence relations allow us to transform the 2 < q < 3 cases into 
incomplete beta functions with positive definite arguments: 

By(a,b) = B(b,a)-B 1 _ y (b,a) (A3) 
B y (a + 1, b) = -1^ [aB y (a, b) - (1 - y) b y a ] (A4) 

where B(a,b) = B n (a,b) is the complete beta function. For y < (a + l)/(a + 6 + 2), B y (a,b) is evaluated using rapidly 
converging continued fractions. For y > (o + l)/(a + i» + 2), equation (A3) is used to transform the problem back into 
a regime where the continued fractions converge quickly. 

For the specific case that b = (e.g. NFW profiles), the transformations cannot be used. However, the continued 
fraction method yields results accurate to better than 1CP 6 for all interesting values of y so long as a is limited to 
physically plausible values (i.e. is limited to 1 to 3, corresponding n = to 2). 

B. CALCULATION OF THE PSF COEFFICIENTS 

Here we discuss the approximations used for calculating the scattering of light from an annulus with inner and outer 
radii (Rj_\,Rj) into an annulus with inner and outer radii (Ri-\,Ri). The general expression for Cv(x, y), the counts 
observed using a detector with a monochromatic point spread function p v (x,y,x',y') observing a monochromatic 
source with flux distribution I v {x, y) is 

C v {x,y)= [ [ I„(x',y')p v (x,y,x',y')dx'dy' (Bl) 



Approximating the source and the PSF as azimuthally symmetric, the above expression becomes 

C V {R)= / R'I v {R')p v {R,R 2 + R' 2 -2RR' cos <p)dR'd(f>, (B2) 
Jo Ja 

where R 2 + R' 2 — 2RR' cos <p is the square of the distance between (x, y) and (x', y') in polar coordinates. Splitting up 
the source flux I(R') into N homogeneous annuli, we have 

C„(R) = y2l v (R j -i,R j ) / R'l{R')p v (R,R 2 + R' 2 -2RR' cos 4>)dR'd4> (B3) 

j=l Jo JRj-i 

where R = (Ri-i + Ri)/2. Each set of N double integrals yields the PSF correction coefficients at energy hi/. We 
evaluate the double integral for four photon energies (hu = 0.37,0.75,1.5,3.6 keV), and interpolate to obtain the 
correction at arbitrary energy. 

C. TECHNICAL DETAILS OF THE REDUCTION AND ANALYSIS 

JACO vl.O includes scripts which prepare the raw data from the Chandra and XMM-Newton archives for analysis. 
These preprocessing scripts yield a uniformly reduced set of spectra for an entire cluster with minimal human inter- 
action. The use of these scripts is not required for undertaking the JACO analysis. Observers may undertake any 
reduction procedure in preparation for JACO, subject to the following requirements on the final output spectra: 

1. Each spectrum must be extracted from a circular annulus surrounding the cluster center. Future versions of 
JACO will also handle elliptical annuli. 

2. The spectrum keyword BACKSCAL must be set to the net source area in square arcminutes (i.e. excluding the 
area of any excised regions intersecting the annulus). 

3. The redistribution matrix file (RMF) and ancilliary response file (ARF) must have the same energy binning for 
all spectra. The largest bin size that preserves relevant calibration features should be used. 

C.l. Pipeline reprocessing 

To begin, the JACO preprocessing scripts reapply the latest version of the standard pipeline analysis to the raw data. 
For Chandra, this has the effect of correcting for the spatial gradients in the filter contamination. For ACIS front- 
illuminated chips with the appropriate detector temperature, we include the charge transfer inefficiency correction 
(CTI) (Townslcy et al. 2000). For XMM-Newton, we rerun the standard pipeline to include the necessary files for 
correcting the EPIC pn spectra for the out-of-time (OOT) events. 
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C.2. Point source removal 

Because we are interested only in the diffuse X-ray emitting medium, we remove contaminating point sources from 
the data. The CIAO wavelet source detection algorithm yields source lists for both the Chandra and the XMM-Newton 
total band images. Visual inspection and correction of the source lists is necessary to remove spurious detections along 
chip borders. For each detected point source, we remove an elliptical region with major axis equal to three times 
the PSF width from photon event files. The source detection program specifies the orientation and axis ratio of the 
ellipses. 

C.3. Flaring background 

Flares in the X-ray background can contaminate the spectra and must therefore be screened using the field total 
light curve at energies > 10 kcV. In many XMM-Newton observations, a substantial fraction of the telescope live time 
is affected by flares. To screen the XMM-Newton data, we apply 2<r clipping of the light curve binned in 100s intervals. 
Then we apply the technique of Nevalainen, Markevitch, & Lumb (2005). They find that small flares significantly affect 
the cr-clipped XMM-Newton data, and that these flares are best filtered by accepting only periods with count rates 
less than 120% of the mean cr-clipped count rate. For Chandra data, we follow the flare rejection method provided as 
a standard contributed background analysis package (Markevitch 2005). 

C.4. Quiescent particle background 

Two types of background affect Chandra and XMM-Newton observations of diffuse sources. The non X-ray particle 
background is dominated by the interaction of charged and non-charged particles with the CCD. In both observatories, 
the particle background comprises a continuum as well as spectral line features whose location depends on the particular 
instrument. Being non X-ray in origin, the particle background is not vignetted — i.e., it is not affected by the variations 
in the effective area across the CCD. To subtract this background from the spectra, we use the latest public versions 
of the blank sky observations for XMM-Newton (Read & Ponman 2003) and Chandra (Markevitch 2005). The 
subtraction is made possible by the fact that after flare rejection, there is little change in the spectral shape of the 
particle background with position or time; only the normalization varies appreciably. 

To carry out the subtraction, we subject the XMM-Newton blank sky spectra to same flare rejection algorithm 
as the cluster spectra (the Chandra blank sky are already calibrated in the same manner as the observations). We 
normalize the spectra of the blank sky fields to match the spectra of the cluster fields observed in the 10-12 keV 
(XMM-Newton) and 8-10 keV (Chandra) energy range (where the effective area for X-ray photons is low). For the 
purposes of calculating the renormalization factor, we exclude emission from the central 100 kpc of the cluster. For 
each extracted cluster spectrum, we use the blank sky fields to extract a matching particle background spectrum. 

C.5. Residual X-ray Background 

Once point sources are removed, the remaining diffuse emission consists of two components: unresolved extragalactic 
AGN (with index ~ 1.4 power law spectra) and unabsorbed, thermal ~ 0.2 keV background with approximately solar 
abundance (Markevitch et al. 2003). Unlike the particle background, the soft X-ray background varies with position 
on the sky, so that it cannot be corrected using blank sky fields. 

One option for dealing with the residual background is to ignore it completely by cutting all emission below 1 keV. 
However, this also removes a significant portion of the cluster signal, including potentially important low-energy line 
emission. Another option is to measure the spectrum in a source-free region and subtract it from the cluster (Read 
& Ponman 2003). However clusters of galaxies often fill the entire field of view of the X-ray detector; it is difficult to 
find a truly source- free region. 

Our approach is to model the residual background as part of the fitting process (e.g. Mahdavi et al. 2005). This 
is the only realistic option when the cluster fills the entire field of view. For each instrument, we add a different 
unabsorbed thermal plasma with free temperature and abundance and arbitrary (negative or positive) normalization. 
To account for the diffuse extragalactic background, we add an absorbed power law with slope fixed to 1.4 and similarly 
free normalization. We require the temperature and metal abundance to be the same for all instruments. However, 
different normalizations are required for each detector because the instrumental background from the blank sky fields 
is different for each instrument. When the renormalized instrumental background is subtracted, the net (positive or 
negative) cosmic background in each field will be different, even though the true cosmic background is the same. 

C.6. Spectrum Extraction 

At present the JACO method includes only spherically symmetric analysis; we therefore extract spectra in concentric 
annuli around the cluster center. The cluster center is determined using the peak of the Chandra X-ray emission, and 
verified using Hubble Space Telescope archival data. 

JACO requires a minimum annulus width of 5". At these widths, the PSF correction (§2.6) becomes important 
even for Chandra data. The widths are by default set to be the same for all the instruments considered. The JACO 
user can set the minimum number of background-subtracted counts required per annulus. The background measured 
at this stage (and this stage only) is roughly estimated from the region outside the largest circular annulus that will 
completely fit in the instruments' field of view. The spectra are grouped by default at 40 events per bin, but the user 
can modify this. 
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To model the spectra, an accurate knowledge of the total solid angle subtended by each annulus is necessary. 
Because an annulus can intersect chip gaps, bad pixels, and excised sources, the calculation of this area is not trivial. 
In particular, the XMM-Newton SAS backs cale procedure cannot calculate accurate areas for the annuli with the 
smallest (5") width. For this reason, JACO creates its own high resolution bad pixel map for each of the three EPIC 
instruments, and uses it to determine the correct solid angles. 

We use the appropriate CIAO and SAS tasks to create response matrices for the spectra. For Chandra, we use 
mkwarf /mkacisrmf ; for XMM-Newton, we use arf gen/rmgf en. The response functions are emission-weighted through 
the use of weight maps. Arnaud et al. (2001) argue that this technique is imperfect, because the spatial variation 
of the uncalibrated source can yield biases in the X-ray temperatures. Their solution, evigweight, is an alternative 
technique for the reduction of XMM-Newton spectra. In evigweight the individual photon events are weighted by 
the inverse of the effective area at the photons' recorded positions, and only one central ARF is created for the entire 
cluster. Because CIAO does not offer a similar procedure for Chandra spectra, we have compared XMM-Newton 
spectra derived using the standard method with spectra derived using evigweight. We find that for the fine binning 
used by JACO, the difference between the standard and evigweight spectra are negligible as long as the detector map 
is extracted for energies where the raw counts spectrum is closest to constant, specifically the channels corresponding 
to the 0.8-1.4 keV photon energy range. 

C.7. Covariant Error Analysis: the Hrothgar Parallel Minimizer 

The command-line (as opposed to the Sherpa-based) frontend to JACO allows the global x 2 to be minimized within 
a clustered computing environment. In essence, Hrothgar runs the Levenberg-Marquard (LM) minimization algorithm 
(Press et al. 1992) as implemented by the GNU Scientific Library (GSL; http://www.gnu.org/software/gsl/). The 
LM algorithm requires knowledge of the Jacobian of the function to be minimized, which is not analytic for X-ray 
spectra due to the highly nonlinear nature of the cooling function. Hrothgar calculates the Jacobian numerically 
using the methods described in Press et al. (1992), with error control as implemented by the GSL library. Because 
calculation of the Jacobian matrix J during the minimization phase requires 2m evaluations of the JACO model for 
m free parameters, each function evaluation can be assigned to a cluster node, resulting in a speed increase directly 
proportional to the number of CPUs, up to 2m CPUs. 

We carry out the error analysis by calculating the covariance matrix C = (JJ T ) _1 . When computing J numerically 
for this purpose, we use both the 3-point and the 5-point numerical differentiation rules. The difference between the 3- 
and the 5-point rules yields an estimate of the error, which we minimize as a function of the step size. The covariance 
matrix obtained in this manner records the interdependence of all the physical parameters in the JACO analysis. The 
confidence contours arc plotted by the methods described in Press et al. (1992). 

Hrothgar will soon be freely downloadable at http://hrothgar.sourceforge.net. 
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